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SIMULATION  OF  CONFINED  PRIMITIVE  ELECTROLYTES:  APPLICATION 
OF  A  NEW  METHOD  OF  SUMMING  THE  COULOMB  FIELD 


Lianrui  Zhang,  Henry  S.  White  and  H.  Ted  Davis 

Department  of  Chemical  Engineering  and  Materials  Science 
University  of  Minnesota 
Minneapolis,  MN  55455 

Abstract 

Recently,  Lekner  has  presented  a  new  method  to  sum  the  Coulomb  forces  between 
charged  particles  of  a  central  system  and  its  images  extended  periodically  in  2  and  3- 
dimensions.  In  this  paper  we  apply  the  new  method  in  canonical  ensemble  Monte  Carlo 
(CMC)  simulations  of  the  primitive  electrolyte  confined  between  two  planar  surfaces;  one 
is  charged  and  the  other  is  neutral.  The  anions  and  cations  have  identical  size  with 
diameter  d  =  4.25  A  and  interact  with  a  hard  sphere  repulsion  and  Coulomb  interaction. 
In  Lekner’s  method  the  long  range  Coulomb  potential  is  computed  from  a  series  of  Bessel 
functions.  We  have  demonstrated  that  the  series  converges  after  about  10  terms  and  so  is 
computationally  simpler  than  the  Ewald  sum  method.  In  our  simulations,  we  obtained  the 
density  distributions  and  mean  electrostatic  potentials  of  the  confined  system  for  the  1:1 
electrolyte  having  concentrations  equal  to  those  of  IM  and  2M  bulk  electrolyte  and  having 
different  surface  charge  densities.  For  large  separation  of  confining  walls,  the  canonical 
ensemble  Monte  Caxlo  results  agree  with  previously  reported  grand  canonical  Monte  Carlo 
results. 
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I.  Introduction 


The  investigation  of  electrolytes  in  confined  geometry  has  attracted  a  lot  of  attention 
in  recent  years.  The  work  has  been  driven  by  the  potential  for  applications  in  electrochem¬ 
istry,  microelectronics,  biological  systems,  colloidal  dispersions,  and  clay  and  soil  systems. 
However,  these  kinds  of  systems  are  difficult  to  study  at  the  molecular  level,  either  ex¬ 
perimentally  or  theoretically.  Therefore  computer  simulation,  usually  either  molecular 
dynamics  or  the  Monte  Carlo  method,  is  a  valuable  tool  for  trying  to  understand  ion 
distributions  and  the  double  layer  potential  of  confined  electrolytes.  Compared  with  the 
minimum  image  method,  the  accuracy  of  a  simulation  is  increased  by  periodic  extension 
of  the  system  in  its  unconfined  dimensions. 

The  common  examples  of  electrical  double  layers  are  electrolytes  near  metal  surfaces 
or  surrounding  the  particles  of  an  electrically  stabilized  colloidal  suspension,  or  the  charge 
distribution  around  the  biological  membranes.  For  point  charges  near  a  planar  surface, 
Gouy  and  Chapman  developed  a  formal  theory  70  years  ago  [l]  to  describe  the  charge 
distributions  near  the  surface  by  using  the  Poisson- Boltzmann  equation  (PBE).  With  the 
realization  of  the  importance  of  double  layers  in  stabilizing  the  colloidal  systems,  the  DLVO 
theory  [2]  aroused  a  renaissance  of  double  layer  theory.  Because  of  the  advent  of  large  fast 
computers  and  the  advancement  of  density  functional  theories,  much  research  has  been 
done  on  double  layers  in  recent  years  [3,4]. 

For  meaningful  computer  simulations,  one  must  have  good  models  of  the  interparti¬ 
cle  interactions.  In  the  case  of  nonpolar  fluids,  interactions  are  short-ranged  and  so  the 
introduction  of  cutoffs  can  be  used  to  reduce  the  cost  of  computations.  However,  the 
coulomb  interactions  of  electrolytes  are  long-ranged  and  so  it  is  not  desirable  to  use  a  cut¬ 
off.  Although  for  bulk  electrolyte  systems  the  minimum  image  short  cut  has  provided  an 
approximate  method  to  account  for  the  long-range  forces,  it  may  cause  unphysical  results 
for  confined  systems.  An  alternative  is  to  use  the  conventional  Ewald  sum  method  [5],  but 
rather  lengthy  summations  must  be  evaluated  by  this  method. 

Recently,  Lekner  has  presented  a  summation  method  that  evaluates  the  long-range 
Coulomb  interactions  in  bulk  or  confined  electrolytes.  The  interaction  potential  is  given  as 
a  series  expansion  of  modified  Bessel  functions  /\;,[6,7].  The  series  is  expected  to  require 
evaluation  of  far  fewer  terms  than  does  the  Ewald  method.  To  examine  the  applicability 
of  Lekner ’s  method  we  have  carried  out  canonical  ensemble  Monte  Carlo  (CMC)  sim\i- 
lations  on  confined  1:1  electrolytes  having  an  average  charge  density  equal  to  bulk  ion 
concentrations  of  IM  and  2M,  respectively.  Different  surface  charge  densities  and  different 
separations  of  the  confining  walls  are  studied.  The  results  are  compared  to  those  of  a 
grand  canonical  Monte  Carlo  (GCMC)  simulation  in  which  the  long-ranged  interactions 
are  handled  by  the  minimum  image  and  the  contributions  from  the  charge  distributions 
everywhere  outside  the  ion’s  minimum  image  [3].  To  further  test  the  scheme  on  confined 
systems,  we  also  carried  out  an  MD  simulation  to  compare  with  the  Ewald  sum  method 
reported  by  Halley  et  al.  [8].  At  the  same  accuracy,  we  found  Leckner's  method  runs 


faster. 


II.  Summation  of  Long  Range  Potentials 

We  give  here  a  brief  sketch  of  Lekner's  summation  method  [7].  Consider  a  simulation 
box  with  side  length  L  along  the  x  and  y  directions  and  H  along  the  r  direction  with  N 
particles  inside.  To  imitate  the  real  system  that  is  confined  between  two  planar  surfaces, 
wc  assume  the  box  is  repeating  to  infinity  in  both  x  and  y  axis.  The  coulomb  interaction 
potential  in  the  medium  with  a  dielectric  constant  e  is  then 


all  boxes 

KJ 


QiQj 

e|rj  —  r 


(1) 


Here  the  summations  over  i  etnd  j  are  from  1  to  N.  If  we  define  the  relative  distances  in 
dimensionless  quantities  (^,77,  and  Cn 


^  =  {xi-  Xj)/L,  T]  =  (y,  -  yj)/L,  (^  =  (z,  -  Zj)IL,  (2) 

with  ^  <  1,T?  <  1,  and  Q  <  H/L 
Then  the  potential  can  be  expressed  as 


OO 

C(r„)=  Y. 


g.gj _ \ _ 

eL  [(^  +  /)2  +  (^  +  rn)2  +C"p/2- 


(3) 


To  sum  up  the  potential  one  uses  the  definition  of  the  F  function 


1 


Tiu) 


oc 

/ 


(d  ) 


the  identity 

OO  I -  OO 

^exp{-{^  +  lft}  =  Wy  ^exp(-7r2/Vt)cos(27r/0,  (o) 

“OO  ^  “OO 

and  the  integral  representations  of  the  modified  Bessel  function  K„ 

^  I 

dt  exp(  — —  rn^t)  =  2(7r|  —  l)*'A'^(27r|/m|).  (6) 

1  ^ 

After  lengthy  mathematical  manipulations,  the  potential  on  each  particle  in  the  central 
box  is  given  in  units  of  q,qj/eL.  by  the  expression 

OO  OO 

q,C)  =  4  cos(27r,^)  y^  A'o(27rm[(r;  +  m  -f  C^])  -  log(coHh(2rrr7)  -  cos( -2C))-  (  ~ » 

/=!  -OO 
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The  modified  Bessel  function  A'o(x)  decays  very  feist  with  distance  [6].  For  an  accuracy  of 
order  10"“*  we  only  need  terms  up  to  x  =  8.  If  we  keep  terms  up  to  x  =  10  the  contributions 
from  the  long  range  tail  will  be  of  order  ~  10“®.  This  is  acceptable  in  the  light  of  the 
simulation  requirement  which  has  been  discussed  by  others  [8,9]. 


III.  Canonical  Ensemble  Monte  Carlo  Simulation 


In  this  paper  we  used  the  canonical  ensemble  Monte  Carlo  to  carry  out  the  simulations. 
The  long  range  coulomb  potential  was  calculated  by  the  method  give  in  Section  II  and  the 
calculations  was  accurate  to  higher  order  term  must  be  <  10“"*.  The  simulation  was  done 
at  different  surface  charge  densities  and  concentrations  for  a  1:1  electrolyte.  Both  ion 
species  are  hard  spheres  with  diameter  d  =  4.25  A.  The  surface  at  z  =  0  is  charged  and 
the  one  at  2  =  if  is  neutral.  The  temperature  for  the  simulation  was  fixed  at  300°K.  The 
parameters  used  in  this  paper  are  given  in  Table  I  for  reference.  In  the  table  AN  is  the 
difference  of  number  of  ions  inside  the  box.  and  N-  are  the  number  of  positive  and 
negative  ions  with  AN  —  N-  —  N+  =  a  A,  where  A  =  L^.  no  is  the  bulk  density  of  the 
electrolyte  which  is  the  asymptotic  value  of  the  electrolyte  and  fixed  to  be  1  M  or  2  M  and 
V  is  the  volume  of  the  box.  In  this  paper,  no  =  0.046/d^  or  0.092/d^. 

During  this  simulation,  the  system  was  equilibrated  for  20,000  steps  to  reach  equilib¬ 
rium  and  another  2  x  10®  steps  to  accumulate  configurations  to  calculate  average  quantities. 
The  density  profiles  emd  mean  electrostatic  potential  0  were  calculated.  At  each  step,  the 
particle  is  moved  to  a  random  position  inside  a  sphere  of  radius  R.  The  move  is  accepted 
if  the  energy  change  is  negative,  i.e.,  SV„m  <  0,  and  accepted  with  the  probability 

exp(-j3Vn)  .  , 

p(n-*m)  =  - —  — r  =  exp(-^Slnm)-  (S) 

exp(-/9Vm) 

when  SV„^  >  0,  where  n  and  m  are  initial  and  final  states  in  an  attempt  and  3  =  1/ kT  with 
k  the  Boltzmann  factor  amd  T  the  temperature.  To  accept  a  move  with  this  probability, 
a  random  number  9.  is  generated  in  the  range  of  0  to  1.  If  ^  <  p(n  — ►  m)  the  move  is 
accepted.  To  sample  the  phase  space  as  much  as  possible  we  vary  R  so  that  the  probability 
of  acceptance  is  about  50  percent. 

Every  10th  of  the  configuration  was  saved  during  the  run  to  accumulate  the  averages 
for  the  calculation  of  density  profiles  n{z)  and  mean  electrostatic  potential  U’(z)  where  c 
is  the  distance  from  the  charged  wall.  The  density  n{z)  was  obtained  by  counting  the 
number  of  particles  in  slices  parallel  to  the  charged  wail  and  dividing  by  the  volume  of  tlie 
slice,  then  is  calculated  via.  the  relation 

47r  /■" 

0(2)=  _/  dziiz  -  Zi)^q,n,{zi).  |9. 
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IV.  Results  and  Discussion 


The  computed  density  profiles  and  mean  electrostatic  potential  are  reported  in  what 
follows.  Surface  charge  density  and  mean  electrostatic  potential  are  reported  in  the  di¬ 
mensionless  forms 


<7 


* 


and  ii>*  =  3ev. 


(10) 


Figure  1  shows  the  density  profiles  from  the  CMC  simulation  compared  with  the  results 
from  the  GCMC  simulations  for  the  surface  charge  density  a*  —  0.42.  The  GCMC  calcu¬ 
lations  are  for  an  open  system  in  equilibrium  with  a  bulk  phase  at  1  M  concentration.  The 
figure  shows  the  density  in  reduced  unit  n,{z)/n,Q,  with  i  specifying  the  ion  species  and 
riiO  is  the  bulk  density  corresponding  to  a  concentration  of  IM  .  From  the  figure  we  see 
that  the  results  of  CMC  and  GCMC  are  in  agreement  within  the  variance  of  the  simulated 
results.  The  mean  electrostatic  potential  from  our  simulations  is  shown  in  figure  2.  We 
have  no  GCMC  results  to  compare  with  the  electrostatic  potential. 


Figures  3  and  4  are  the  density  and  mean  electrostatic  potential  profiles  for  the  case 
of  a*  =  0.7  and  H  =  30d.  Except  for  the  small  shift  of  the  second  peak  in  the  density 
profile,  CMC  and  GCMC  results  are  quite  close  and  both  show  the  layer  of  negative  ions 
at  ^  ~  3d/2.  The  mean  electrostatic  potentials  of  the  CMC  and  GCMC  calculations  are 
also  in  good  agreement. 


To  see  how  the  density  profiles  and  the  mean  electrostatic  potentials  would  vary  at 
higher  concentration  A  CMC  simulation  was  carried  out  at  a  2M  concentration  of  a  1:1 
electrolyte  with  a*  =  0.396  and  H  =  I3d.  These  results  are  shown  in  Figures  5  and  6.  We 
can  see  that  the  ion  densities  approach  the  bulk  values  faster  than  in  the  1  M  concentration 
electrolyte.  The  ion  densities  predicted  by  the  CMC  and  GCMC  simulations  agree  well, 
but  the  ion  density  profiles  differ  sufficiently  to  yield  observable  differences  in  the  mean 
electrostatic  potential.  The  diffuse  layer  potential  xi'*(^d)  and  the  total  potential  drop 
b'*(0)  from  CMC  and  GCMC  eure  given  in  Table  II  for  comparison.  The  errors  given 
in  Table  II  are  obtained  from  the  difference  between  the  mean  value  of  half  simulation 
averages  with  the  final  average  results.  We  do  not  know  if  the  small  differences  arise 
from  real  differences  between  ensembles  or  from  differences  between  Lekner's  sum  and  the 
minimum  image  potential  used  in  the  GCMC  simulations. 
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Table  I.  Parcimeters  used  in  this  paper,  a*  =  acP  je.  side  length  £,  surface  separation 
number  of  positive  and  negative  ions  and  iV_,  —  Nj^  =  a  A,  A  =  £-,  no 

is  the  bulk  concentration  in  unit  l/d^. 


Hid) 

Lid) 

(7* 

AV 

N- 

AA' 

no( l/d^) 

29.0 

3.09 

0.42 

12 

16 

4 

0.046 

30.0 

2.93 

0.7 

11 

17 

6 

0.046 

13.0 

4.50 

0.396 

21 

29 

8 

0.092 

Table  11.  The  diffuse  layer  potential  and  the  total  potential  drop  v^(0)  for 

different  surface  charge  densities  a*  and  bulk  concentration  r?o  from  CMC  and  GCMC. 


riz) 

a* 

Uo 

CMC 

GCMC 

ri}d) 

0.42 

0.046 

3.13(0.13) 

3.08(0.10) 

0.7 

0.046 

5.28(0.04) 

5.71(0.14) 

0.396 

0.092 

1.83(0.006) 

2.29(0.09) 

4’*i0) 

0.42 

0.046 

7.58(0.13) 

7.52 

0.7 

0.046 

12.70(0.04) 

13.10 

0.396 

0.092 

6.15(0.006) 

6.47 
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Figure  Captions 


Figure  1.  Density  profiles  n[z)/no  for  the  1;1  electrolyte  at  surface  charge  density  of  a*  — 
ad}  jt  =  0.42  and  separation  H  ~  29c/  and  concentration  1  M.  Circles  and  black  dots 
for  GCMC,  and  triangles  are  for  CMC. 

Figure  2.  Mean  electrostatic  potential  w*(z)  for  the  1:1  electrolyte  for  the  CMC  simulation  of 
Figure  1. 

Figure  3.  Density  profiles  n(z)/no  for  the  same  conditions  as  in  figure  1  except  charge  density 
of  a*  =  0.7  and  separation  H  ~  30c/.  Solid  lines  for  MGC  results,  circles  and  black 
dots  are  for  MCMC,  and  triangles  are  for  CMC. 

Figure  4.  Mean  electrostatic  potential  rp*{z)  for  the  1:1  electrolyte  at  the  same  conditions  as  in 
figure  3.  The  dotted  line  corresponds  to  CMC  results,  and  the  dashed  lines  to  GCMC 
results. 

Figure  5.  Density  profiles  n{z)/no  for  the  1:1  electrolyte  at  surface  charge  density  of  a*  —  0.396 
and  separation  H  ~  13c/  and  concentration  2  M.  Circles  and  black  dots  denote  MCMC 
results,  and  triangles  denote  CMC  results. 

Figure  6.  Mean  electrostatic  potential  ip*{z)  for  the  1:1  electrolyte  at  the  same  conditions  as 
in  figure  5.  Dotted  line  correspond  to  CMC  results,  and  dashed  lines  correspond  to 
GCMC  results. 
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